home *** CD-ROM | disk | FTP | other *** search
/ STraTOS 1997 April & May / STraTOS 1 - 1997 April & May.iso / CD01 / INTERNET / SITES / LITTLE / P3SRC.ZIP / ATARI / SOR.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-04-25  |  25.0 KB  |  1,227 lines
  1. /****************************************************************************
  2. *                   sor.c
  3. *
  4. *  This module implements functions that manipulate surfaces of revolution.
  5. *
  6. *  This module was written by Dieter Bayer [DB].
  7. *
  8. *  from Persistence of Vision(tm) Ray Tracer
  9. *  Copyright 1996 Persistence of Vision Team
  10. *---------------------------------------------------------------------------
  11. *  NOTICE: This source code file is provided so that users may experiment
  12. *  with enhancements to POV-Ray and to port the software to platforms other
  13. *  than those supported by the POV-Ray Team.  There are strict rules under
  14. *  which you are permitted to use this file.  The rules are in the file
  15. *  named POVLEGAL.DOC which should be distributed with this file. If
  16. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  17. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  18. *  Forum.  The latest version of POV-Ray may be found there as well.
  19. *
  20. * This program is based on the popular DKB raytracer version 2.12.
  21. * DKBTrace was originally written by David K. Buck.
  22. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  23. *
  24. *****************************************************************************/
  25.  
  26. /****************************************************************************
  27. *
  28. *  Explanation:
  29. *
  30. *    The surface of revolution primitive is defined by a set of points
  31. *    in 2d space wich are interpolated by cubic splines. The resulting
  32. *    2d function is rotated about an axis to form the final object.
  33. *
  34. *    All calculations are done in the object's (u,v,w)-coordinate system
  35. *    with the (w)-axis being the rotation axis.
  36. *
  37. *    One spline segment in the (r,w)-plane is given by the equation
  38. *
  39. *      r^2 = f(w) = A * w * w * w + B * w * w + C * w + D.
  40. *
  41. *    To intersect a ray R = P + k * D transformed into the object's
  42. *    coordinate system with the surface of revolution, the equation
  43. *
  44. *      (Pu + k * Du)^2 + (Pv + k * Dv)^2 = f(Pw + k * Dw)
  45. *
  46. *    has to be solved for k (cubic polynomial in k).
  47. *
  48. *    Note that Pu, Pv, Pw and Du, Dv, Dw denote the coordinates
  49. *    of the vectors P and D.
  50. *
  51. *  Syntax:
  52. *
  53. *    revolution
  54. *    {
  55. *      number_of_points,
  56. *
  57. *      <P[0]>, <P[1]>, ..., <P[n-1]>
  58. *
  59. *      [ open ]
  60. *    }
  61. *
  62. *    Note that the P[i] are 2d vectors where u corresponds to the radius
  63. *    and v to the height.
  64. *
  65. *    Note that the first and last point, i.e. P[0] and P[n-1], are used
  66. *    to determine the derivatives at the end point.
  67. *
  68. *    Note that the x coordinate of a point corresponds to the radius and
  69. *    the y coordinate to the height; the z coordinate isn't used.
  70. *
  71. *  ---
  72. *
  73. *  Ideas for the surface of revolution were taken from:
  74. *
  75. *    P. Burger and D. Gillies, "Rapid Ray Tracing of General Surfaces
  76. *    of Revolution", New Advances in Computer Graphics, Proceedings
  77. *    of CG International '89, R. A. Earnshaw, B. Wyvill (Eds.),
  78. *    Springer, ..., pp. 523-531
  79. *
  80. *  ---
  81. *
  82. *  May 1994 : Creation. [DB]
  83. *
  84. *****************************************************************************/
  85.  
  86. #include "frame.h"
  87. #include "povray.h"
  88. #include "vector.h"
  89. #include "povproto.h"
  90. #include "bbox.h"
  91. #include "polysolv.h"
  92. #include "matrices.h"
  93. #include "objects.h"
  94. #include "sor.h"
  95.  
  96.  
  97.  
  98. /*****************************************************************************
  99. * Local preprocessor defines
  100. ******************************************************************************/
  101.  
  102. /* Minimal intersection depth for a valid intersection. */
  103.  
  104. #define DEPTH_TOLERANCE 1.0e-4
  105.  
  106. /* |Coefficients| < COEFF_LIMIT are regarded to be 0. */
  107.  
  108. #define COEFF_LIMIT 1.0e-20
  109.  
  110. /* Part of the surface of revolution hit. */
  111.  
  112. #define BASE_PLANE 1
  113. #define CAP_PLANE  2
  114. #define CURVE      3
  115.  
  116.  
  117.  
  118. /*****************************************************************************
  119. * Local typedefs
  120. ******************************************************************************/
  121.  
  122. typedef struct Sor_Intersection_Structure SOR_INT;
  123.  
  124. struct Sor_Intersection_Structure
  125. {
  126.   DBL d;
  127.   int t, n;
  128. };
  129.  
  130.  
  131.  
  132. /*****************************************************************************
  133. * Static functions
  134. ******************************************************************************/
  135.  
  136. static int intersect_sor PARAMS((RAY *Ray, SOR *Sor, SOR_INT *Int));
  137. static int  All_Sor_Intersections PARAMS((OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack));
  138. static int  Inside_Sor PARAMS((VECTOR point, OBJECT *Object));
  139. static void Sor_Normal PARAMS((VECTOR Result, OBJECT *Object, INTERSECTION *Inter));
  140. static void *Copy_Sor PARAMS((OBJECT *Object));
  141. static void Translate_Sor PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  142. static void Rotate_Sor PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  143. static void Scale_Sor PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  144. static void Transform_Sor PARAMS((OBJECT *Object, TRANSFORM *Trans));
  145. static void Invert_Sor PARAMS((OBJECT *Object));
  146. static void Destroy_Sor PARAMS((OBJECT *Object));
  147.  
  148.  
  149. /*****************************************************************************
  150. * Local variables
  151. ******************************************************************************/
  152.  
  153. static METHODS Sor_Methods =
  154. {
  155.   All_Sor_Intersections,
  156.   Inside_Sor, Sor_Normal,
  157.   Copy_Sor,
  158.   Translate_Sor, Rotate_Sor,
  159.   Scale_Sor, Transform_Sor, Invert_Sor, Destroy_Sor
  160. };
  161.  
  162.  
  163.  
  164. /*****************************************************************************
  165. *
  166. * FUNCTION
  167. *
  168. *   All_Sor_Intersections
  169. *
  170. * INPUT
  171. *
  172. *   Object      - Object
  173. *   Ray         - Ray
  174. *   Depth_Stack - Intersection stack
  175. *   
  176. * OUTPUT
  177. *
  178. *   Depth_Stack
  179. *   
  180. * RETURNS
  181. *
  182. *   int - TRUE, if a intersection was found
  183. *   
  184. * AUTHOR
  185. *
  186. *   Dieter Bayer
  187. *   
  188. * DESCRIPTION
  189. *
  190. *   Determine ray/surface of revolution intersection and
  191. *   clip intersection found.
  192. *
  193. * CHANGES
  194. *
  195. *   May 1994 : Creation.
  196. *
  197. ******************************************************************************/
  198.  
  199. static int All_Sor_Intersections(Object, Ray, Depth_Stack)
  200. OBJECT *Object;
  201. RAY *Ray;
  202. ISTACK *Depth_Stack;
  203. {
  204.   int i, max_i, Found;
  205.   SOR_INT *I;
  206.   VECTOR IPoint;
  207.  
  208.   /* Allocate memory for intersections. */
  209.  
  210.   I = (SOR_INT *)POV_MALLOC(5*((SOR *)Object)->Number*sizeof(SOR_INT), "surface of revolution intersection array");
  211.  
  212.   Found = FALSE;
  213.  
  214.   max_i = intersect_sor(Ray, (SOR *)Object, I);
  215.  
  216.   if (max_i)
  217.   {
  218.     for (i = 0; i < max_i; i++)
  219.     {
  220.       VEvaluateRay(IPoint, Ray->Initial, I[i].d, Ray->Direction);
  221.  
  222.       if (Point_In_Clip(IPoint, Object->Clip))
  223.       {
  224.         push_entry_i1_i2(I[i].d,IPoint,Object,I[i].t,I[i].n,Depth_Stack);
  225.  
  226.         Found = TRUE;
  227.       }
  228.     }
  229.   }
  230.  
  231.   POV_FREE (I);
  232.  
  233.   return(Found);
  234. }
  235.  
  236.  
  237.  
  238. /*****************************************************************************
  239. *
  240. * FUNCTION
  241. *
  242. *   intersect_sor
  243. *
  244. * INPUT
  245. *
  246. *   Ray          - Ray
  247. *   Sor   - Sor
  248. *   Intersection - Sor intersection structure
  249. *   
  250. * OUTPUT
  251. *
  252. *   Intersection
  253. *   
  254. * RETURNS
  255. *
  256. *   int - Number of intersections found
  257. *   
  258. * AUTHOR
  259. *
  260. *   Dieter Bayer
  261. *   
  262. * DESCRIPTION
  263. *
  264. *   Determine ray/surface of revolution intersection.
  265. *
  266. *   NOTE: The curve is rotated about the y-axis!
  267. *         Order reduction cannot be used.
  268. *
  269. * CHANGES
  270. *
  271. *   May 1994 : Creation.
  272. *
  273. ******************************************************************************/
  274.  
  275. static int intersect_sor(Ray, Sor, Intersection)
  276. RAY *Ray;
  277. SOR *Sor;
  278. SOR_INT *Intersection;
  279. {
  280.   int i = 0, j, n;
  281.   DBL a, b, k, l, len, u, v, x[4], y[3], r0;
  282.   VECTOR P, D;
  283.   SOR_SPLINE_ENTRY Entry;
  284.  
  285.   Increase_Counter(stats[Ray_Sor_Tests]);
  286.  
  287.   /* Init intersection structure. */
  288.  
  289.   Intersection->d = 0.0;
  290.   Intersection->n =
  291.   Intersection->t = 0;
  292.  
  293.   /* Transform the ray into the surface of revolution space. */
  294.  
  295.   MInvTransPoint(P, Ray->Initial, Sor->Trans);
  296.  
  297.   MInvTransDirection(D, Ray->Direction, Sor->Trans);
  298.  
  299.   VLength(len, D);
  300.  
  301.   VInverseScaleEq(D, len);
  302.  
  303.   /* Test if ray misses object's bounds. */
  304.  
  305. #ifdef SOR_EXTRA_STATS
  306.   Increase_Counter(stats[Sor_Bound_Tests]);
  307. #endif
  308.  
  309.   if (((D[Y] >= 0.0) && (P[Y] >  Sor->Height2)) ||
  310.       ((D[Y] <= 0.0) && (P[Y] <  Sor->Height1)) ||
  311.       ((D[X] >= 0.0) && (P[X] >  Sor->Radius2)) ||
  312.       ((D[X] <= 0.0) && (P[X] < -Sor->Radius2)))
  313.   {
  314.     return(FALSE);
  315.   }
  316.  
  317.   /* Get distance of the ray from rotation axis (= y axis). */
  318.  
  319.   r0 = P[X] * D[Z] - P[Z] * D[X];
  320.  
  321.   if ((a = D[X] * D[X] + D[Z] * D[Z]) > 0.0)
  322.   {
  323.     r0 /= sqrt(a);
  324.   }
  325.  
  326.   /* Test if ray misses object's bounds. */
  327.  
  328.   if (r0 > Sor->Radius2)
  329.   {
  330.     return(FALSE);
  331.   }
  332.  
  333. #ifdef SOR_EXTRA_STATS
  334.   Increase_Counter(stats[Sor_Bound_Tests_Succeeded]);
  335. #endif
  336.  
  337.   if (fabs(D[Y]) > DEPTH_TOLERANCE)
  338.   {
  339.     /* Test base/cap plane. */
  340.  
  341.     if (Test_Flag(Sor, CLOSED_FLAG))
  342.     {
  343.       /* Test base plane. */
  344.  
  345.       if (Sor->Base_Radius_Squared > DEPTH_TOLERANCE)
  346.       {
  347.         a = (Sor->Height1 - P[Y]) / D[Y];
  348.  
  349.         if ((a > DEPTH_TOLERANCE) && (a < Max_Distance))
  350.         {
  351.           u = P[X] + a * D[X];
  352.           v = P[Z] + a * D[Z];
  353.  
  354.           b = u * u + v * v;
  355.  
  356.           if (b <= Sor->Base_Radius_Squared)
  357.           {
  358.             Intersection[i].t   = BASE_PLANE;
  359.             Intersection[i++].d = a / len;
  360.           }
  361.         }
  362.       }
  363.  
  364.       /* Test cap plane. */
  365.  
  366.       if (Sor->Cap_Radius_Squared > DEPTH_TOLERANCE)
  367.       {
  368.         a = (Sor->Height1 - P[Y]) / D[Y];
  369.  
  370.         if ((a > DEPTH_TOLERANCE) && (a < Max_Distance))
  371.         {
  372.           u = P[X] + a * D[X];
  373.           v = P[Z] + a * D[Z];
  374.  
  375.           b = u * u + v * v;
  376.  
  377.           if (b <= Sor->Cap_Radius_Squared)
  378.           {
  379.             Intersection[i].t   = BASE_PLANE;
  380.             Intersection[i++].d = a / len;
  381.           }
  382.         }
  383.       }
  384.     }
  385.   }
  386.  
  387.   /* Test all segments. */
  388.  
  389.   for (j = 0; j < Sor->Number; j++)
  390.   {
  391.     Entry = Sor->Spline->Entry[j];
  392.  
  393.     /* Test if ray misses segment's bounds. */
  394.  
  395. #ifdef SOR_EXTRA_STATS
  396.     Increase_Counter(stats[Sor_Bound_Tests]);
  397. #endif
  398.  
  399.     if (((D[Y] >= 0.0) && (P[Y] >  Entry.h2)) ||
  400.         ((D[Y] <= 0.0) && (P[Y] <  Entry.h1)) ||
  401.         ((D[X] >= 0.0) && (P[X] >  Entry.r2)) ||
  402.         ((D[X] <= 0.0) && (P[X] < -Entry.r2)))
  403.     {
  404.       continue;
  405.     }
  406.  
  407.     if (r0 > Entry.r2)
  408.     {
  409.       continue;
  410.     }
  411.  
  412. #ifdef SOR_EXTRA_STATS
  413.     Increase_Counter(stats[Sor_Bound_Tests_Succeeded]);
  414. #endif
  415.  
  416.     /* Cubic curve. */
  417.  
  418.     x[0] = Entry.A * D[Y] * D[Y] * D[Y];
  419.  
  420.     x[1] = D[Y] * D[Y] * (3.0 * Entry.A * P[Y] + Entry.B) - D[X] * D[X] - D[Z] * D[Z];
  421.  
  422.     x[2] = D[Y] * (P[Y] * (3.0 * Entry.A * P[Y] + 2.0 * Entry.B) + Entry.C) - 2.0 * (P[X] * D[X] + P[Z] * D[Z]);
  423.  
  424.     x[3] = P[Y] * (P[Y] * (Entry.A * P[Y] + Entry.B) + Entry.C) + Entry.D - P[X] * P[X] - P[Z] * P[Z];
  425.  
  426.     n = Solve_Polynomial(3, x, y, Test_Flag(Sor, STURM_FLAG), 0.0);
  427.  
  428.     while (n--)
  429.     {
  430.       k = y[n];
  431.  
  432.       if ((k > DEPTH_TOLERANCE) && (k < Max_Distance))
  433.       {
  434.         l = P[Y] + k * D[Y];
  435.  
  436.         if ((l >= Entry.h1) && (l <= Entry.h2))
  437.         {
  438.           Intersection[i].t   = CURVE;
  439.           Intersection[i].n   = j;
  440.           Intersection[i++].d = k / len;
  441.         }
  442.       }
  443.     }
  444.   }
  445.  
  446.   if (i)
  447.   {
  448.     Increase_Counter(stats[Ray_Sor_Tests_Succeeded]);
  449.   }
  450.  
  451.   return(i);
  452. }
  453.  
  454.  
  455.  
  456. /*****************************************************************************
  457. *
  458. * FUNCTION
  459. *
  460. *   Inside_Sor
  461. *
  462. * INPUT
  463. *
  464. *   IPoint - Intersection point
  465. *   Object - Object
  466. *   
  467. * OUTPUT
  468. *   
  469. * RETURNS
  470. *
  471. *   int - TRUE if inside
  472. *   
  473. * AUTHOR
  474. *
  475. *   Dieter Bayer
  476. *   
  477. * DESCRIPTION
  478. *
  479. *   Return true if point lies inside the surface of revolution.
  480. *
  481. * CHANGES
  482. *
  483. *   May 1994 : Creation.
  484. *
  485. ******************************************************************************/
  486.  
  487. static int Inside_Sor(IPoint, Object)
  488. VECTOR IPoint;
  489. OBJECT *Object;
  490. {
  491.   int i;
  492.   DBL r0, r;
  493.   VECTOR P;
  494.   SOR *Sor = (SOR *)Object;
  495.   SOR_SPLINE_ENTRY Entry;
  496.  
  497.   /* Transform the point into the surface of revolution space. */
  498.  
  499.   MInvTransPoint(P, IPoint, Sor->Trans);
  500.  
  501.   /* Test if we are inside the cylindrical bound. */
  502.  
  503.   if ((P[Y] >= Sor->Height1) && (P[Y] <= Sor->Height2))
  504.   {
  505.     r0 = P[X] * P[X] + P[Z] * P[Z];
  506.  
  507.     /* Test if we are inside the cylindrical bound. */
  508.  
  509.     if (r0 <= Sqr(Sor->Radius2))
  510.     {
  511.       /* Now find the segment the point is in. */
  512.  
  513.       for (i = 0; i < Sor->Number; i++)
  514.       {
  515.         Entry = Sor->Spline->Entry[i];
  516.  
  517.         /* P[Y] >= R->Entry[i].h1 needn't be tested. */
  518.  
  519.         if ((P[Y] >= Entry.h1) && (P[Y] <= Entry.h2))
  520.         {
  521.           break;
  522.         }
  523.       }
  524.  
  525.       /* Have we found any segment? */
  526.  
  527.       if (i < Sor->Number)
  528.       {
  529.         r = P[Y] * (P[Y] * (P[Y] * Entry.A + Entry.B) + Entry.C) + Entry.D;
  530.  
  531.         if (r0 <= r)
  532.         {
  533.           /* We're inside. */
  534.  
  535.           return(!Test_Flag(Sor, INVERTED_FLAG));
  536.         }
  537.       }
  538.     }
  539.   }
  540.  
  541.   /* We're outside. */
  542.  
  543.   return(Test_Flag(Sor, INVERTED_FLAG));
  544. }
  545.  
  546.  
  547.  
  548. /*****************************************************************************
  549. *
  550. * FUNCTION
  551. *
  552. *   Sor_Normal
  553. *
  554. * INPUT
  555. *
  556. *   Result - Normal vector
  557. *   Object - Object
  558. *   Inter  - Intersection found
  559. *   
  560. * OUTPUT
  561. *
  562. *   Result
  563. *
  564. * RETURNS
  565. *   
  566. * AUTHOR
  567. *
  568. *   Dieter Bayer
  569. *   
  570. * DESCRIPTION
  571. *
  572. *   Calculate the normal of the surface of revolution in a given point.
  573. *
  574. * CHANGES
  575. *
  576. *   May 1994 : Creation.
  577. *
  578. ******************************************************************************/
  579.  
  580. static void Sor_Normal(Result, Object, Inter)
  581. OBJECT *Object;
  582. VECTOR Result;
  583. INTERSECTION *Inter;
  584. {
  585.   DBL k;
  586.   VECTOR P;
  587.   SOR *Sor = (SOR *)Object;
  588.   SOR_SPLINE_ENTRY Entry;
  589.   VECTOR N;
  590.  
  591.   Make_Vector(N, 0.0, 1.0, 0.0);
  592.  
  593.   if (Inter->i1 == CURVE)
  594.   {
  595.     /* Transform the intersection point into the surface of revolution space. */
  596.  
  597.     MInvTransPoint(P, Inter->IPoint, Sor->Trans);
  598.  
  599.     if (P[X] * P[X] + P[Z] * P[Z] > DEPTH_TOLERANCE)
  600.     {
  601.       Entry = Sor->Spline->Entry[Inter->i2];
  602.  
  603.       k = 0.5 * (P[Y] * (3.0 * Entry.A * P[Y] + 2.0 * Entry.B) + Entry.C);
  604.  
  605.       N[X] = P[X];
  606.       N[Y] = -k;
  607.       N[Z] = P[Z];
  608.     }
  609.   }
  610.  
  611.   /* Transform the normalt out of the surface of revolution space. */
  612.  
  613.   MTransNormal(Result, N, Sor->Trans);
  614.  
  615.   VNormalize(Result, Result);
  616. }
  617.  
  618.  
  619.  
  620. /*****************************************************************************
  621. *
  622. * FUNCTION
  623. *
  624. *   Translate_Sor
  625. *
  626. * INPUT
  627. *
  628. *   Object - Object
  629. *   Vector - Translation vector
  630. *   
  631. * OUTPUT
  632. *
  633. *   Object
  634. *   
  635. * RETURNS
  636. *   
  637. * AUTHOR
  638. *
  639. *   Dieter Bayer
  640. *   
  641. * DESCRIPTION
  642. *
  643. *   Translate a surface of revolution.
  644. *
  645. * CHANGES
  646. *
  647. *   May 1994 : Creation.
  648. *
  649. ******************************************************************************/
  650.  
  651. static void Translate_Sor(Object, Vector, Trans)
  652. OBJECT *Object;
  653. VECTOR Vector;
  654. TRANSFORM *Trans;
  655. {
  656.   Transform_Sor(Object, Trans);
  657. }
  658.  
  659.  
  660.  
  661. /*****************************************************************************
  662. *
  663. * FUNCTION
  664. *
  665. *   Rotate_Sor
  666. *
  667. * INPUT
  668. *
  669. *   Object - Object
  670. *   Vector - Rotation vector
  671. *   
  672. * OUTPUT
  673. *
  674. *   Object
  675. *   
  676. * RETURNS
  677. *   
  678. * AUTHOR
  679. *
  680. *   Dieter Bayer
  681. *   
  682. * DESCRIPTION
  683. *
  684. *   Rotate a surface of revolution.
  685. *
  686. * CHANGES
  687. *
  688. *   May 1994 : Creation.
  689. *
  690. ******************************************************************************/
  691.  
  692. static void Rotate_Sor(Object, Vector, Trans)
  693. OBJECT *Object;
  694. VECTOR Vector;
  695. TRANSFORM *Trans;
  696. {
  697.   Transform_Sor(Object, Trans);
  698. }
  699.  
  700.  
  701.  
  702. /*****************************************************************************
  703. *
  704. * FUNCTION
  705. *
  706. *   Scale_Sor
  707. *
  708. * INPUT
  709. *
  710. *   Object - Object
  711. *   Vector - Scaling vector
  712. *   
  713. * OUTPUT
  714. *
  715. *   Object
  716. *   
  717. * RETURNS
  718. *   
  719. * AUTHOR
  720. *
  721. *   Dieter Bayer
  722. *   
  723. * DESCRIPTION
  724. *
  725. *   Scale a surface of revolution.
  726. *
  727. * CHANGES
  728. *
  729. *   May 1994 : Creation.
  730. *
  731. ******************************************************************************/
  732.  
  733. static void Scale_Sor(Object, Vector, Trans)
  734. OBJECT *Object;
  735. VECTOR Vector;
  736. TRANSFORM *Trans;
  737. {
  738.   Transform_Sor(Object, Trans);
  739. }
  740.  
  741.  
  742.  
  743. /*****************************************************************************
  744. *
  745. * FUNCTION
  746. *
  747. *   Transform_Sor
  748. *
  749. * INPUT
  750. *
  751. *   Object - Object
  752. *   Trans  - Transformation to apply
  753. *   
  754. * OUTPUT
  755. *
  756. *   Object
  757. *   
  758. * RETURNS
  759. *   
  760. * AUTHOR
  761. *
  762. *   Dieter Bayer
  763. *   
  764. * DESCRIPTION
  765. *
  766. *   Transform a surface of revolution and recalculate its bounding box.
  767. *
  768. * CHANGES
  769. *
  770. *   May 1994 : Creation.
  771. *
  772. ******************************************************************************/
  773.  
  774. static void Transform_Sor(Object, Trans)
  775. OBJECT *Object;
  776. TRANSFORM *Trans;
  777. {
  778.   Compose_Transforms(((SOR *)Object)->Trans, Trans);
  779.  
  780.   Compute_Sor_BBox((SOR *)Object);
  781. }
  782.  
  783.  
  784.  
  785. /*****************************************************************************
  786. *
  787. * FUNCTION
  788. *
  789. *   Invert_Sor
  790. *
  791. * INPUT
  792. *
  793. *   Object - Object
  794. *   
  795. * OUTPUT
  796. *
  797. *   Object
  798. *   
  799. * RETURNS
  800. *   
  801. * AUTHOR
  802. *
  803. *   Dieter Bayer
  804. *   
  805. * DESCRIPTION
  806. *
  807. *   Invert a surface of revolution.
  808. *
  809. * CHANGES
  810. *
  811. *   May 1994 : Creation.
  812. *
  813. ******************************************************************************/
  814.  
  815. static void Invert_Sor(Object)
  816. OBJECT *Object;
  817. {
  818.   Invert_Flag(Object, INVERTED_FLAG);
  819. }
  820.  
  821.  
  822.  
  823. /*****************************************************************************
  824. *
  825. * FUNCTION
  826. *
  827. *   Create_Sor
  828. *
  829. * INPUT
  830. *   
  831. * OUTPUT
  832. *   
  833. * RETURNS
  834. *
  835. *   SOR * - new surface of revolution
  836. *   
  837. * AUTHOR
  838. *
  839. *   Dieter Bayer
  840. *   
  841. * DESCRIPTION
  842. *
  843. *   Create a new surface of revolution.
  844. *
  845. * CHANGES
  846. *
  847. *   May 1994 : Creation.
  848. *
  849. ******************************************************************************/
  850.  
  851. SOR *Create_Sor()
  852. {
  853.   SOR *New;
  854.  
  855.   New = (SOR *)POV_MALLOC(sizeof(SOR), "surface of revolution");
  856.  
  857.   INIT_OBJECT_FIELDS(New,SOR_OBJECT,&Sor_Methods)
  858.  
  859.   New->Trans = Create_Transform();
  860.  
  861.   New->Spline = NULL;
  862.  
  863.   New->Radius2             =
  864.   New->Base_Radius_Squared =
  865.   New->Cap_Radius_Squared  = 0.0;
  866.  
  867.   return(New);
  868. }
  869.  
  870.  
  871.  
  872. /*****************************************************************************
  873. *
  874. * FUNCTION
  875. *
  876. *   Copy_Sor
  877. *
  878. * INPUT
  879. *
  880. *   Object - Object
  881. *   
  882. * OUTPUT
  883. *   
  884. * RETURNS
  885. *
  886. *   void * - New surface of revolution
  887. *   
  888. * AUTHOR
  889. *
  890. *   Dieter Bayer
  891. *   
  892. * DESCRIPTION
  893. *
  894. *   Copy a surface of revolution structure.
  895. *
  896. *   NOTE: The splines are not copied, only the number of references is
  897. *         counted, so that Destray_Sor() knows if they can be destroyed.
  898. *
  899. * CHANGES
  900. *
  901. *   May 1994 : Creation.
  902. *
  903. *   Sep 1994 : fixed memory leakage [DB]
  904. *
  905. ******************************************************************************/
  906.  
  907. static void *Copy_Sor(Object)
  908. OBJECT *Object;
  909. {
  910.   SOR *New, *Sor = (SOR *)Object;
  911.  
  912.   New = Create_Sor();
  913.  
  914.   /* Get rid of the transformation created in Create_Sor(). */
  915.  
  916.   Destroy_Transform(New->Trans);
  917.  
  918.   /* Copy surface of revolution. */
  919.  
  920.   *New = *Sor;
  921.  
  922.   New->Trans = Copy_Transform(Sor->Trans);
  923.  
  924.   New->Spline->References++;
  925.  
  926.   return(New);
  927. }
  928.  
  929.  
  930.  
  931. /*****************************************************************************
  932. *
  933. * FUNCTION
  934. *
  935. *   Destroy_Sor
  936. *
  937. * INPUT
  938. *
  939. *   Object - Object
  940. *   
  941. * OUTPUT
  942. *
  943. *   Object
  944. *   
  945. * RETURNS
  946. *   
  947. * AUTHOR
  948. *
  949. *   Dieter Bayer
  950. *   
  951. * DESCRIPTION
  952. *
  953. *   Destroy a surface of revolution.
  954. *
  955. *   NOTE: The splines are destroyed if they are no longer used by any copy.
  956. *
  957. * CHANGES
  958. *
  959. *   May 1994 : Creation.
  960. *
  961. ******************************************************************************/
  962.  
  963. static void Destroy_Sor (Object)
  964. OBJECT *Object;
  965. {
  966.   SOR *Sor = (SOR *)Object;
  967.  
  968.   Destroy_Transform(Sor->Trans);
  969.  
  970.   if (--(Sor->Spline->References) == 0)
  971.   {
  972.     POV_FREE (Sor->Spline->Entry);
  973.  
  974.     POV_FREE (Sor->Spline);
  975.   }
  976.  
  977.   POV_FREE (Object);
  978. }
  979.  
  980.  
  981.  
  982. /*****************************************************************************
  983. *
  984. * FUNCTION
  985. *
  986. *   Compute_Sor_BBox
  987. *
  988. * INPUT
  989. *
  990. *   Sor - Sor
  991. *   
  992. * OUTPUT
  993. *
  994. *   Sor
  995. *   
  996. * RETURNS
  997. *   
  998. * AUTHOR
  999. *
  1000. *   Dieter Bayer
  1001. *   
  1002. * DESCRIPTION
  1003. *
  1004. *   Calculate the bounding box of a surface of revolution.
  1005. *
  1006. * CHANGES
  1007. *
  1008. *   May 1994 : Creation.
  1009. *
  1010. ******************************************************************************/
  1011.  
  1012. void Compute_Sor_BBox(Sor)
  1013. SOR *Sor;
  1014. {
  1015.   Make_BBox(Sor->BBox, -Sor->Radius2, Sor->Height1, -Sor->Radius2,
  1016.     2.0 * Sor->Radius2, Sor->Height2 - Sor->Height1, 2.0 * Sor->Radius2);
  1017.  
  1018.   Recompute_BBox(&Sor->BBox, Sor->Trans);
  1019. }
  1020.  
  1021.  
  1022.  
  1023. /*****************************************************************************
  1024. *
  1025. * FUNCTION
  1026. *
  1027. *   Compute_Sor
  1028. *
  1029. * INPUT
  1030. *
  1031. *   Sor - Sor
  1032. *   P          - Points defining surface of revolution
  1033. *   
  1034. * OUTPUT
  1035. *
  1036. *   Sor
  1037. *   
  1038. * RETURNS
  1039. *   
  1040. * AUTHOR
  1041. *
  1042. *   Dieter Bayer, June 1994
  1043. *   
  1044. * DESCRIPTION
  1045. *
  1046. *   Calculate the spline segments of a surface of revolution
  1047. *   from a set of points.
  1048. *
  1049. *   Note that the number of points in the surface of revolution has to be set.
  1050. *
  1051. * CHANGES
  1052. *
  1053. *   May 1994 : Creation.
  1054. *
  1055. ******************************************************************************/
  1056.  
  1057. void Compute_Sor(Sor, P)
  1058. SOR *Sor;
  1059. UV_VECT *P;
  1060. {
  1061.   int i, n;
  1062.   DBL A, B, C, D, w, k[4];
  1063.   DBL x[4], xmax;
  1064.   DBL y[2], ymax, ymin;
  1065.   DBL c[3], r[2];
  1066.   MATRIX Mat;
  1067.  
  1068.   /* Allocate Sor->Number segments. */
  1069.  
  1070.   if (Sor->Spline == NULL)
  1071.   {
  1072.     Sor->Spline = (SOR_SPLINE *)POV_MALLOC(sizeof(SOR_SPLINE), "spline segments of surface of revoluion");
  1073.  
  1074.       Sor->Spline->References = 1;
  1075.  
  1076.     Sor->Spline->Entry = (SOR_SPLINE_ENTRY *)POV_MALLOC(Sor->Number*sizeof(SOR_SPLINE_ENTRY), "spline segments of surface of revoluion");
  1077.   }
  1078.   else
  1079.   {
  1080.     Error("Surface of revolution segments are already defined.\n");
  1081.   }
  1082.  
  1083.   /* We want to know the size of the overall bounding cylinder. */
  1084.  
  1085.   xmax = ymax = -BOUND_HUGE;
  1086.          ymin =  BOUND_HUGE;
  1087.  
  1088.   /* Calculate segments, i.e. cubic patches. */
  1089.  
  1090.   for (i = 0; i < Sor->Number; i++)
  1091.   {
  1092.     if ((fabs(P[i+2][Y] - P[i][Y]) < EPSILON) ||
  1093.         (fabs(P[i+3][Y] - P[i+1][Y]) < EPSILON))
  1094.     {
  1095.       Error("Incorrect point in surface of revolution.\n");
  1096.     }
  1097.  
  1098.     /* Use cubic interpolation. */
  1099.  
  1100.     k[0] = P[i+1][X] * P[i+1][X];
  1101.     k[1] = P[i+2][X] * P[i+2][X];
  1102.     k[2] = (P[i+2][X] - P[i][X]) / (P[i+2][Y] - P[i][Y]);
  1103.     k[3] = (P[i+3][X] - P[i+1][X]) / (P[i+3][Y] - P[i+1][Y]);
  1104.  
  1105.     k[2] *= 2.0 * P[i+1][X];
  1106.     k[3] *= 2.0 * P[i+2][X];
  1107.  
  1108.     w = P[i+1][Y];
  1109.  
  1110.     Mat[0][0] = w * w * w;
  1111.     Mat[0][1] = w * w;
  1112.     Mat[0][2] = w;
  1113.     Mat[0][3] = 1.0;
  1114.  
  1115.     Mat[2][0] = 3.0 * w * w;
  1116.     Mat[2][1] = 2.0 * w;
  1117.     Mat[2][2] = 1.0;
  1118.     Mat[2][3] = 0.0;
  1119.  
  1120.     w = P[i+2][Y];
  1121.  
  1122.     Mat[1][0] = w * w * w;
  1123.     Mat[1][1] = w * w;
  1124.     Mat[1][2] = w;
  1125.     Mat[1][3] = 1.0;
  1126.  
  1127.     Mat[3][0] = 3.0 * w * w;
  1128.     Mat[3][1] = 2.0 * w;
  1129.     Mat[3][2] = 1.0;
  1130.     Mat[3][3] = 0.0;
  1131.  
  1132.     MInvers(Mat, Mat);
  1133.  
  1134.     /* Calculate coefficients of cubic patch. */
  1135.  
  1136.     A = k[0] * Mat[0][0] + k[1] * Mat[0][1] + k[2] * Mat[0][2] + k[3] * Mat[0][3];
  1137.     B = k[0] * Mat[1][0] + k[1] * Mat[1][1] + k[2] * Mat[1][2] + k[3] * Mat[1][3];
  1138.     C = k[0] * Mat[2][0] + k[1] * Mat[2][1] + k[2] * Mat[2][2] + k[3] * Mat[2][3];
  1139.     D = k[0] * Mat[3][0] + k[1] * Mat[3][1] + k[2] * Mat[3][2] + k[3] * Mat[3][3];
  1140.  
  1141.     if (fabs(A) < EPSILON) A = 0.0;
  1142.     if (fabs(B) < EPSILON) B = 0.0;
  1143.     if (fabs(C) < EPSILON) C = 0.0;
  1144.     if (fabs(D) < EPSILON) D = 0.0;
  1145.  
  1146.     Sor->Spline->Entry[i].A = A;
  1147.     Sor->Spline->Entry[i].B = B;
  1148.     Sor->Spline->Entry[i].C = C;
  1149.     Sor->Spline->Entry[i].D = D;
  1150.  
  1151.     /* Get maximum radius**2 in current segment. */
  1152.  
  1153.     y[0] = P[i+1][Y];
  1154.     y[1] = P[i+2][Y];
  1155.  
  1156.     x[0] = x[2] = P[i+1][X];
  1157.     x[1] = x[3] = P[i+2][X];
  1158.  
  1159.     c[0] = 3.0 * A;
  1160.     c[1] = 2.0 * B;
  1161.     c[2] = C;
  1162.  
  1163.     n = Solve_Polynomial(2, c, r, FALSE, 0.0);
  1164.  
  1165.     while (n--)
  1166.     {
  1167.       if ((r[n] >= y[0]) && (r[n] <= y[1]))
  1168.       {
  1169.         x[n] = sqrt(r[n] * (r[n] * (r[n] * A + B) + C) + D);
  1170.       }
  1171.     }
  1172.  
  1173.     /* Set current segment's bounding cylinder. */
  1174.  
  1175.     Sor->Spline->Entry[i].h1 = y[0];
  1176.     Sor->Spline->Entry[i].h2 = y[1];
  1177.  
  1178.     Sor->Spline->Entry[i].r2 = max(max(x[0], x[1]), max(x[2], x[3]));
  1179.  
  1180.     /* Keep track of overall bounding cylinder. */
  1181.  
  1182.     xmax = max(xmax, Sor->Spline->Entry[i].r2);
  1183.  
  1184.     ymin = min(ymin, Sor->Spline->Entry[i].h1);
  1185.     ymax = max(ymax, Sor->Spline->Entry[i].h2);
  1186.   }
  1187.  
  1188.   /* Set overall bounding cylinder. */
  1189.  
  1190.   Sor->Radius2 = xmax;
  1191.  
  1192.   Sor->Height1 = ymin;
  1193.   Sor->Height2 = ymax;
  1194.  
  1195.   /* Get cap radius. */
  1196.  
  1197.   w = Sor->Spline->Entry[Sor->Number-1].h1;
  1198.  
  1199.   A = Sor->Spline->Entry[Sor->Number-1].A;
  1200.   B = Sor->Spline->Entry[Sor->Number-1].B;
  1201.   C = Sor->Spline->Entry[Sor->Number-1].C;
  1202.   D = Sor->Spline->Entry[Sor->Number-1].D;
  1203.  
  1204.   if ((Sor->Cap_Radius_Squared = w * (w * (A * w + B) + C) + D) < 0.0)
  1205.   {
  1206.     Sor->Cap_Radius_Squared = 0.0;
  1207.   }
  1208.  
  1209.   /* Get base radius. */
  1210.  
  1211.   w = Sor->Spline->Entry[0].h1;
  1212.  
  1213.   A = Sor->Spline->Entry[0].A;
  1214.   B = Sor->Spline->Entry[0].B;
  1215.   C = Sor->Spline->Entry[0].C;
  1216.   D = Sor->Spline->Entry[0].D;
  1217.  
  1218.   if ((Sor->Base_Radius_Squared = w * (w * (A * w + B) + C) + D) < 0.0)
  1219.   {
  1220.     Sor->Base_Radius_Squared = 0.0;
  1221.   }
  1222. }
  1223.  
  1224.  
  1225.  
  1226.  
  1227.